1.1.1 Lancement de R studio

Sauf à être complètement masochiste, on n’utilise jamais R directement mais on lance d’abord l’interface R-Studio qui facilite conisdérablement l’ensemble des opérations et offre une gamme considérable de services. Il ne faut toutefois pas confondre les deux et il serait par exemple ridicule d’indiquer sur un CV en vue d’un emploi de statisticien que l’on sait utiliser R-studio en oubliant de préciser que l’on maîtrise R.

1.1.2 Création d’un projet

Pour créer un projet on utilise le menus déroulant File/new project/ … et on définit un dossier de notre ordinateur (existant ou à créer) qui contiendra le projet. Une fois l’opération effectuée, on pourra constater que ce dossier contient un fichier xxx.Rproj ou xxx est en principe le nom du dossier dans lequel vous avez stocké le projet.

Ce fichier contient toute une série d’informations dont nous ne parlerons pas dans ce cours d’initiation mais qui, pour faire simple, définissent un ensemble de réglages du logiciel et de préférences de l’utilisateur.

Si vous fermez Rstudio (faites-le !) il vous suffira pour reprendre votre travail là où vous vous étiez arrêté :

• de lancer R-studio et de cliquer sur File/open project/… suivi du nom du fichier xxx.Rproj
• ou plus simplement encore de double cliquer sur le fichier xxx.Rproj ce qui lancera automatiquement Rstudio

Le dossier contenant votre projet R peut être organisé à votre convenance. Certains mettent tout les fichier pêle-mêle dans le dossier. D’autres préfèrent créer des sous-dossiers contenant des données, des programmes, des résultats, des figures. Vous déciderez à l’usage de ce qui vous convient le mieux, mais le point important est que tout ce qui entre ou sort de vos programmes R doit être de préférence stocké dans le répertoire du projet.

knitr::include_graphics("figures/monprojet.png")

1.2.1 Round 1. Saisie des données et affichage du tableau

On crée un programme R avec File/New File/R Script puis on l’enregistre avec File/Save/ … suivi du nom du programme.

# Saisie des variables
  CODE <- c("A","B","C","D","E")
  X <- c(10,20,40,50,180)
  Y <- c(40,60,40,60,50)
  
  # Regroupement dans un tableau
  coo <- data.frame(X,Y)
  
  # Ajout du nom des lignes
  row.names(coo) <- CODE
  
  # Affichage du tableau
  coo

    X  Y
  A  10 40
  B  20 60
  C  40 40
  D  50 60
  E 180 50

Normalement, les étudiants qui utilisent un tableur ont du aller plus vite et Excel mène sur R par 1-0

1.2.2 Round 2. Affichage de la carte

Vous devez essayez de reproduire l’image correspondant au problème posé

plot(X,Y, 
       col="red", 
       pch=20, 
       xlim=c(0,180),
       ylim=c(0,90),
       asp = 1)
  text(X,Y,
       labels = CODE, 
       pos = 2)

La création d’un graphique est à première vue plus facile avec un logiciel click-bouton. L’avantage est très clairement pour Excel qui mène désormais 2 à 0.

1.2.3 Round 3. Calcul de la station la plus accessible à vol d’oiseau (distance euclidienne)

Vous devez calculer une matrice de distance euclidienne entre toutes les stations et trouver la plus accessible.

# calcul la matrice de distance euclidienne
  mat<-dist(coo, upper = T, method = "euclidean")
  mat

          A         B         C         D         E
  A            22.36068  30.00000  44.72136 170.29386
  B  22.36068            28.28427  30.00000 160.31220
  C  30.00000  28.28427            22.36068 140.35669
  D  44.72136  30.00000  22.36068           130.38405
  E 170.29386 160.31220 140.35669 130.38405          

# distance moyenne
  apply(as.matrix(mat),1,mean)

        A         B         C         D         E 
   53.47518  48.19143  44.20033  45.49322 120.26936 

Là, je parie que les utilisateurs d’Excel ont eu un peu plus de mal … En tous les cas, Excel ne mèneplus que par 2 à 1

1.2.4 Round 4. Calcul de la station la plus accessible par la route (distance de Manhattan)

Vous devez calculer une matrice de distance de Manhattan entre toutes les stations et trouver la plus accessible.

# calcul la matrice de distance de Manhattan
  mat<-dist(coo,upper = T, method = "manhattan")
  mat

    A   B   C   D   E
  A      30  30  60 180
  B  30      40  30 170
  C  30  40      30 150
  D  60  30  30     140
  E 180 170 150 140    

# distance moyenne de Manhattan
  apply(as.matrix(mat),1,mean)

  A   B   C   D   E 
   60  54  50  52 128 

Je reconnais que c’est unpeu facile, mais à nouveau R l’emporte ce qui fait désormais match nul 2-2

1.2.5 Round 5. Localisation du dépôt de carburant

Dans le cas particulier de la distance de Manhattan, le calcul du point le plus proche de tous les autres s’obtient facilement en calculant le point médian dont les coordonnées correspondent à la médiane de X et la médiane de Y.

medX <- median(X)
  medX

[1] 40

medY <- median(Y)
  medY

[1] 50

A priori, le calcul est aussi facile dans R et dans Excel : match nul 3-3

1.2.6 Round 6. Localisation de la caserne de pompiers

Dans le cas particulier de la distance de Manhattan, le calcul du point minimisant la distance maximale s’obtient en trouvant le centre du diamètre minimal en X et en Y. Il s’agit de la localisation la plus équitable où le plus défavorisé est le moins défavorisé possible.

equX <- (max(X)+min(X))/2
  equX

[1] 95

equY <- (max(Y)+min(Y))/2
  equY

[1] 50

A priori, le calcul est toujours aussi facile dans R et dans Excel : match nul 4-4

1.2.7 Round 7. Visualisation des deux points sur la carte

On va placer en bleu le point médian et en vert le point le plus équitable. Dans le cas de R on peut recopier les lignes de code du graphique du round n°2 ce qui gagne désormais du temps :

# Programme antérieur
  plot(X,Y, 
       col="red", 
       pch=20, 
       xlim=c(0,180),
       ylim=c(0,90),
       asp = 1)
  text(X,Y,
       labels = CODE, 
       pos = 2)
  
  # Ajout du dépôt de carburant
  points(medX, medY, col="blue", pch=3)
  text(medX,medY, "dépot",pos=1)
  
  # Ajout du point médian
  points(equX, equY, col="green", pch=3)
  text(equX,equY, "caserne",pos=1)

Le résultat du match est incertain mais R n’est plus désavantagé puisqu’on peut recycler les lignes de code précédentes pour le graphique de base. Disons 5-5 même s’il y a de fortes chances que R l’emporte.

1.2.8 Dernier round. Refaire toute l’analyse avec une station de plus

Deux stations F(100,20) et G(150,30) ont été ajoutées. Il faut refaire la carte finale. Cela ne pose aucun problème dans R puisqu’il suffit de modifier l’entrée des données et récupérer des bouts de programme

# (1) Saisie des variables
  CODE <- c("A","B","C","D","E","F","G")
  X <- c(10,20,40,50,180,100,150)
  Y <- c(40,60,40,60,50,20,30)
  coo <- data.frame(X,Y)
  row.names(coo) <- CODE
  
  
  # (2) calcul des points centraux
  medX <- median(X)
  medY <- median(Y)
  equX <- (max(X)+min(X))/2
  equY <- (max(Y)+min(Y))/2
  
  
  # (3) Graphique
  plot(X,Y, 
       col="red", 
       pch=20, 
       xlim=c(0,180),
       ylim=c(0,90),
       asp = 1)
  text(X,Y,
       labels = CODE, 
       pos = 2)
  
  # Ajout du dépôt de carburant
  points(medX, medY, col="blue", pch=3)
  text(medX,medY, "Dépôt",pos=1)
  
  # Ajout du point médian
  points(equX, equY, col="green", pch=3)
  text(equX,equY, "Caserne",pos=1)

Excel n’a aucune chance d’aller plus vite et R remporte le match par KO !

2.1.1 Définir un espace de travail

le répertoire de travail dans lequel seront lus ou écrit les fichiers.

Dans le programme ci-dessous, les lignes commençant par “#” sont des commentaires. La seule ligne active est getwd() qui indique la position du répertoire courant. Il est recommandé aux débutants, mais aussi aux programmeurs plus expérimentés de ne pas hésiter à ajouter beaucoup de lignes de commentaires. Cela facilité l’apprentissage et permet de transmettre ses programmes à d’autres.

# ---------------------------
  # (1) ESPACE DE TRAVAIL
  # ---------------------------
  # (1.1) Quel est le répertoire actuel ?
  getwd()

[1] "/Users/claudegrasland1/git/Init_R_2021"

# ---------------------------
  # (1) ESPACE DE TRAVAIL
  # ---------------------------
  # (1.1) Quel est le répertoire actuel ?
  getwd()

[1] "C:/git/Init_R_2021"

On repère ensuite le chemin de l’emplacement du dossier où se trouvent les données et on examine son contenu. Attention, contrairement à windows il faut toujours utiliser des slash (/) et non pas des antislash pour décrire les chemins d’accès des fichiers

# (1.2) Choix du rÈpertoire contenant les donnÈes
  monchemin<-"data/euro1988"
  list.files(monchemin)

[1] "euro1988.rda" "euro1988.txt"

Nous allons maintenant essayer de charger le fichier de données

2.2.1 Importation d’un tableau

# --------------------------------------------------------
  # (2) IMPORTATION ET MISE EN FORME D'UN TABLEAU DE DONNEES
  # --------------------------------------------------------
  # (2.1) Importation d'un fichier .txt
  euro <- read.table("data/euro1988/euro1988.txt", 
                     dec=".",
                     header=TRUE)

Pour savoir si le chargement s’est bien pssé on peut afficher le tableau en tapant son nom

euro 

  PAYS BLOC   PNB  TMI ESP URB NAT MOR FEC JEU VIE SUP POP    X    Y
  1  ALB  Soc   600 43.0  71  34  27   6 3.3  35   5  29 3.1 1445  684
  2  AUT  Cap 10000 10.3  75  55  12  12 1.4  18  14  84 7.6 1121 1031
  3  BEL  Cap  9200  9.7  75  95  12  11 1.5  19  14  31 9.9  734 1187
  4  BUL  Soc  2000 14.5  72  65  13  11 2.0  21  11 111 9.0 1649  839
  5  DAN  Cap 12600  8.4  75  84  11  11 1.5  18  15  43 5.1  887 1500
   [ reached 'max' / getOption("max.print") -- omitted 19 rows ]

Mais on peut aussi se contenter d’afficher les premières lignes (head) ou les dernières lignes (tail) du tableau. Par exemple, pour voir les trois premières et les 5 dernières lignes :

head(euro,3)

  PAYS BLOC   PNB  TMI ESP URB NAT MOR FEC JEU VIE SUP POP    X    Y
  1  ALB  Soc   600 43.0  71  34  27   6 3.3  35   5  29 3.1 1445  684
  2  AUT  Cap 10000 10.3  75  55  12  12 1.4  18  14  84 7.6 1121 1031
  3  BEL  Cap  9200  9.7  75  95  12  11 1.5  19  14  31 9.9  734 1187

tail(euro,5)

   PAYS BLOC   PNB  TMI ESP URB NAT MOR FEC JEU VIE SUP  POP    X    Y
  20  ROY  Cap  8900  9.5  75  91  13  12 1.8  19  15 245 57.1  336 1452
  21  SUE  Cap 13200  5.9  77  83  12  11 1.8  18  18 450  8.4 1061 1931
  22  SUI  Cap 17800  6.8  77  61  12   9 1.5  17  14  41  6.6  882  956
  23  TCH  Soc  3200 13.9  71  74  14  12 2.0  24  11 128 15.6 1219 1157
  24  YOU  Soc  2300 27.1  70  47  15   8 2.1  24   8 256 23.6 1356  855

Naturellement … les instructions présentées ci-dessus ne marcheront pas pour tous les tableaux et vous devez vous attendre à pas mal de difficultés à ce stade. Il existe en effet beaucoup d’options à connaître pour charger les fichiers. Vous pouvez ainsi examiner toutes les options de la procédure read.table en tapant un ? suivi du nom de la procédure. Et vous pouvez faire de même pour les autres procédures d’importation comme read.csv, read.delim, etc…

2.2.2 Choix du nom des lignes

Par défaut, le nom des lignes correspond à l’ordre de lecture du fichier (1..24). Mais il peut être intéressant de le préciser en lui attribuant par exemple la valeur d’une variable servant d’identifiant.

rownames(euro)<-euro$PAYS
  head(euro)

    PAYS BLOC   PNB  TMI ESP URB NAT MOR FEC JEU VIE SUP POP    X    Y
  ALB  ALB  Soc   600 43.0  71  34  27   6 3.3  35   5  29 3.1 1445  684
  AUT  AUT  Cap 10000 10.3  75  55  12  12 1.4  18  14  84 7.6 1121 1031
  BEL  BEL  Cap  9200  9.7  75  95  12  11 1.5  19  14  31 9.9  734 1187
  BUL  BUL  Soc  2000 14.5  72  65  13  11 2.0  21  11 111 9.0 1649  839
  DAN  DAN  Cap 12600  8.4  75  84  11  11 1.5  18  15  43 5.1  887 1500
   [ reached 'max' / getOption("max.print") -- omitted 1 rows ]

2.2.3 Vérification du type des variables

On peut s’éviter beaucoup d’ennuis en contrôlant le type des variables qui a été attribué aux différentes colonnes par R lors de la lecture d’un fichier. Pour cela on commence par regarder quels sont les types de variables du tableau à l’aide de l’instruction str() :

str(euro)

'data.frame':   24 obs. of  15 variables:
   $ PAYS: chr  "ALB" "AUT" "BEL" "BUL" ...
   $ BLOC: chr  "Soc" "Cap" "Cap" "Soc" ...
   $ PNB : int  600 10000 9200 2000 12600 4800 12200 10100 3700 2000 ...
   $ TMI : num  43 10.3 9.7 14.5 8.4 9 5.8 8 12.3 19 ...
   $ ESP : int  71 75 75 72 75 76 74 75 74 70 ...
   $ URB : int  34 55 95 65 84 91 62 73 58 58 ...
   $ NAT : int  27 12 12 13 11 12 12 14 11 12 ...
   $ MOR : int  6 12 11 11 11 8 10 10 9 14 ...
   $ FEC : num  3.3 1.4 1.5 2 1.5 1.7 1.6 1.8 1.7 1.8 ...
   $ JEU : int  35 18 19 21 18 23 19 21 21 21 ...
   $ VIE : int  5 14 14 11 15 12 13 13 13 13 ...
   $ SUP : int  29 84 31 111 43 505 337 551 132 93 ...
   $ POP : num  3.1 7.6 9.9 9 5.1 39 4.9 55.9 10.1 10.6 ...
   $ X   : int  1445 1121 734 1649 887 358 1299 646 1453 1341 ...
   $ Y   : int  684 1031 1187 839 1500 582 2103 941 593 1042 ...

L’instruction str() est très importante dans R. Elle permet en effet d’examiner le type d’un objet et des éléments qui le composent. Ici on note que l’objet euro est un data.frame c’est-à-dire un tableau de données croisant des indivdus (lignes) et des variables de types différents (colonnes). Par défaut, R a attribué le type “Factor” à toutes les variables composées de texte, le type “int” ou integer à toutes les variables comportant uniquement des nombres entiers et enfin le type “num” ou numérique à toutes les variables comportant des chiffres avec des décimales. Il peut arriver que ce choix par défaut ne correspondent pas aux traitements que l’on veut effectuer et, dans ce cas, on procèdera à des conversions de variables d’un type à un autre.

Ainsi le type “Factor” est difficile à utiliser pour les débutants et il est souvent confondu avec le type “Character”. Pour fixer les idées, on peut dire qu’un “Factor” correspond à ces classes portant à la fois un numéro et une étiquette. Ainsi la variable “BLOC” est bien un facteur correspondant à deux classes : celle des pays capitalistes (numéro = 1, étiquette =“cap”) et celle des pays socialistes (numero = 2, étiquette = “soc”). Par contre la variable PAYS ne définit pas des classes puisqu’elle est différente pour chaque pays et est donc du type “character”.

De la même manière, un nombre entier peut conduire à des comportements bizarre pour les valeurs élevées et il vaut mieux utiliser le tupe numérique lorsqu’on est en présence de nombre réels qui ont juste été arrondis à zéro chiffres après la virgule. Dans notre exemple, il n’existe pas de variables entières mais uniquement des nombres réels arrondis.

euro$PAYS<-as.character(euro$PAYS)
  euro$BLOC<-as.factor(euro$BLOC)
  euro$PNB<-as.numeric(euro$PNB)
  euro$ESP<-as.numeric(euro$ESP)
  euro$URB<-as.numeric(euro$URB)
  euro$NAT<-as.numeric(euro$NAT)
  euro$MOR<-as.numeric(euro$MOR)
  euro$JEU<-as.numeric(euro$JEU)
  euro$VIE<-as.numeric(euro$VIE)
  euro$SUP<-as.numeric(euro$SUP)
  euro$POP<-as.numeric(euro$POP)
  euro$X<-as.numeric(euro$X)
  euro$Y<-as.numeric(euro$Y)
  str(euro)

'data.frame':   24 obs. of  15 variables:
   $ PAYS: chr  "ALB" "AUT" "BEL" "BUL" ...
   $ BLOC: Factor w/ 2 levels "Cap","Soc": 2 1 1 2 1 1 1 1 1 2 ...
   $ PNB : num  600 10000 9200 2000 12600 4800 12200 10100 3700 2000 ...
   $ TMI : num  43 10.3 9.7 14.5 8.4 9 5.8 8 12.3 19 ...
   $ ESP : num  71 75 75 72 75 76 74 75 74 70 ...
   $ URB : num  34 55 95 65 84 91 62 73 58 58 ...
   $ NAT : num  27 12 12 13 11 12 12 14 11 12 ...
   $ MOR : num  6 12 11 11 11 8 10 10 9 14 ...
   $ FEC : num  3.3 1.4 1.5 2 1.5 1.7 1.6 1.8 1.7 1.8 ...
   $ JEU : num  35 18 19 21 18 23 19 21 21 21 ...
   $ VIE : num  5 14 14 11 15 12 13 13 13 13 ...
   $ SUP : num  29 84 31 111 43 505 337 551 132 93 ...
   $ POP : num  3.1 7.6 9.9 9 5.1 39 4.9 55.9 10.1 10.6 ...
   $ X   : num  1445 1121 734 1649 887 ...
   $ Y   : num  684 1031 1187 839 1500 ...

Et voilà : désormais notre tableau ne comporte plus que des variables de type caractère ou de type numérique. Cela nous a demandé un peu de temps mais nous évitera par la suite beaucoup d’ennuis. En effet, au moins 50 % des erreurs dans R viennent d’un problème de mauvaise définiton du type d’’un objet ou du type d’une variable !!!

2.2.4 Enregistrement du tableau au format propre à R

Il n’est pas forcément necessaire d’enregistrer ses données au format propre à R (.rda)s’il s’agit d’un petit tableau. Il suffit en effet de sauvegarder le programme que nous venons d’écrire et de l’executer chaque fois que nécessaire. Toutefois, il peut arriver que cette étape soit longue et dans ce cas la sauvegarde est utile.

save(euro, file="data/euro1988/euro1988.rda")

on pourra ensuite facilement récupérer les données par la commande suivante

load("data/euro1988/euro1988.rda")

2.2.5 Un petit résumé statistique

N’oublions pas que R est avant tout un logiciel de statistiques… La récompense d’une bonne définition du tableau sera d’obtenir pour chaque variable un résumé adapté à son type :

summary(euro)

     PAYS            BLOC         PNB             TMI              ESP       
   Length:24          Cap:16   Min.   :  600   Min.   : 5.800   Min.   :70.00  
   Class :character   Soc: 8   1st Qu.: 2275   1st Qu.: 8.475   1st Qu.:71.75  
   Mode  :character            Median : 6850   Median : 9.600   Median :74.50  
                               Mean   : 7208   Mean   :13.121   Mean   :73.67  
                               3rd Qu.:10575   3rd Qu.:14.825   3rd Qu.:75.00  
        URB             NAT             MOR             FEC       
   Min.   :30.00   Min.   :10.00   Min.   : 6.00   Min.   :1.400  
   1st Qu.:57.50   1st Qu.:12.00   1st Qu.: 9.00   1st Qu.:1.575  
   Median :68.00   Median :12.50   Median :10.50   Median :1.700  
   Mean   :67.92   Mean   :13.46   Mean   :10.33   Mean   :1.829  
   3rd Qu.:83.25   3rd Qu.:14.00   3rd Qu.:11.00   3rd Qu.:2.000  
        JEU             VIE            SUP             POP       
   Min.   :15.00   Min.   : 5.0   Min.   : 29.0   Min.   : 3.10  
   1st Qu.:19.00   1st Qu.:11.0   1st Qu.: 80.5   1st Qu.: 7.35  
   Median :20.00   Median :13.0   Median :130.0   Median :10.45  
   Mean   :21.33   Mean   :12.5   Mean   :198.5   Mean   :20.64  
   3rd Qu.:23.25   3rd Qu.:14.0   3rd Qu.:304.0   3rd Qu.:27.20  
         X                Y       
   Min.   : 147.0   Min.   : 567  
   1st Qu.: 760.2   1st Qu.: 851  
   Median :1042.0   Median :1100  
   Mean   : 991.7   Mean   :1154  
   3rd Qu.:1309.5   3rd Qu.:1348  
   [ getOption("max.print") est atteint -- ligne 1 omise ]

On va tenter de répondre à des questions simples sans détailler dans l’immédiat les outils statistiques qui permettent de valider scientifiquement les réponses. On insistera plus ici sur les capacités de visualisation de résultats sous R.

2.3.1 Analyse de la répartition des pays par bloc politique (BLOC)

Deux petites commandes pour dénombrer l’effectif des classes et les représenter

# Analyse de la  variable qualitative (BLOC)
  #calcul des fréquences
  summary(euro$BLOC)

Cap Soc 
   16   8 

#diagramme en bâtons
  plot(euro$BLOC)

#diagramme en bâtons
  plot(euro$BLOC)

2.3.2 Analyse de la distribution des taux de mortalité infantile (TMI)

Quelques petites commandes pour camlculer les paramètres principaux de la mortalité infantile et la représenter.

# Analyse de la variable quantitative (TMI)
  
  # calcul et visualisation des quantiles
  summary(euro$TMI)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    5.800   8.475   9.600  13.121  14.825  43.000 

boxplot(euro$TMI)

#calcul de la moyenne et de l'écart type
  mean(euro$TMI)

boxplot(euro$TMI)

#calcul de la moyenne et de l'écart type
  mean(euro$TMI)

[1] 13.12083

sd(euro$TMI)

[1] 8.513875

# visualisation de l'histogramme
  hist(euro$TMI)

# visualisation de l'histogramme
  hist(euro$TMI)

2.3.3 Comaparaison de la mortalité infantile par bloc (BLOC*TMI)

Quelques outils pour analyser conjointement les deux distributions.

# calcul et visualisation des quantiles
  tapply(euro$TMI, euro$BLOC, summary)

$Cap
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    5.800   7.925   8.650   9.069   9.800  15.800 

  $Soc
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
     9.20   14.35   18.25   21.23   25.98   43.00 

boxplot(euro$TMI~euro$BLOC)

#calcul de la moyenne et de l'écart type
  tapply(euro$TMI, euro$BLOC, mean)

boxplot(euro$TMI~euro$BLOC)

#calcul de la moyenne et de l'écart type
  tapply(euro$TMI, euro$BLOC, mean)

     Cap      Soc 
   9.06875 21.22500 

tapply(euro$TMI, euro$BLOC, sd)

      Cap       Soc 
   2.434945 10.624197 

# visualisation des deux histogrammes
  par(mfrow= c(2,1)) 
  hist(euro$TMI[euro$BLOC =="Soc"],
  #     xlim = c(0,50),
       main="Pays Socialistes", 
       col="red") 
  hist(euro$TMI[euro$BLOC =="Cap"],
  #     xlim = c(0,50),
       main="Pays Capitalistes",
       col="blue") 

# visualisation des deux histogrammes
  par(mfrow= c(2,1)) 
  hist(euro$TMI[euro$BLOC =="Soc"],
  #     xlim = c(0,50),
       main="Pays Socialistes", 
       col="red") 
  hist(euro$TMI[euro$BLOC =="Cap"],
  #     xlim = c(0,50),
       main="Pays Capitalistes",
       col="blue") 

2.3.4 Test d’égalité des moyennes (BLOC*TMI)

Maintenant que nous avons pris connaissance de la forme de nos deux variables BLOC et TMI, essayons de réaliser une opération statistique de plus haut niveau : le test d’une hypothèse. Nous savons grâce aux analyses précédentes que la mortalité infantile semble plus élévée en général dans les pays socialistes que dans les pays capitalistes. Mais pouvons nous le *prouver" de façon rigoureuse en utilisant les règles de l’art ?

Bien évidemment, ce n’est pas l’apprentissage de R qui va vous apprendre la théorie des tests statistiques. Mais R va vous permettre très rapidement de mettre en oeuvre les tests que vous aurez appris en cours ou dans les manuels. Et surtout, il va vous offrir une gamme de possibilités de tests et de variantes qui devrait vous pousser à renforcer simultanément vos connaissances en R et en statitiques.

Dans l’exemple qui nous intéresse, les étudiants ayant une formation de base en statistique savent qu’il existe un test d’égalité des moyennes de deux échantillons appelé test t de student qui semble approprié au problème de comapraison de la mortalité infantile des pays socialistes et capitalistes. Il tient compte de l’effectif de ces échantillons (combiens de pays de chaque type ?) et de leur hétérogénéité interne (appelée variance). On devine intuitivement que la différence des moyennes de deux échantillons sera d’autant plus significative qu’elle oppose deux groupes homogènes comportant chacun de nombreux pays.

Si cela ne vous parait pas intuitif, imaginez un premier groupe de trois pays ayant pour mortalité (10, 20, 30) et un second groupe de trois pays également ayant pour mortalité (5,15,25). Leurs moyennes respectives sont 20 et 15, soit un écart de 5. Mais il semblerait hasardeux de conclure à de réelles différences compte tenu de l’étalement de chaque distribution.

# Comparaison de la distribution de deux échantillons
  VARQUANT<-c(10,20,30,5,15,25)
  VARQUALI<-c("A","A","A","B","B","B")
  par(mfrow=c(1,1))
  boxplot(VARQUANT~ VARQUALI)

Imaginez maintenant un premier groupe de 6 pays ayant pour mortalité (18,19,20,20,21,22) et un second groupe de 10 pays ayant pour mortalité (13,14,15,15,16,17). Les moyennes des deux groupes sont toujours respectivement de 20 et de 15, mais on “sent” que les écarts sont plus pertinents car ils portent sur deux groupes plus larges et surtout plus homogènes. Dit autrement, il est peu probable que les différences observées dans la seconde expérience soient le fruit du hasard alors que, dans le premier cas, on pouvait imaginer que les pays des deux groupes n’étaient pas fondamentalement différents.

# Comparaison de la distribution de deux échantillons
  VARQUANT<-c(18,19,20,20,21,22,13,14,15,15,16,17)
  VARQUALI<-c("A","A","A","A","A","A","B","B","B","B","B","B")
  par(mfrow=c(1,1))
  boxplot(VARQUANT~ VARQUALI)

# Comparaison de la distribution de deux échantillons
  VARQUANT<-c(10,20,30,5,15,25)
  VARQUALI<-c("A","A","A","B","B","B")
  par(mfrow=c(1,1))
  boxplot(VARQUANT~ VARQUALI)

Imaginez maintenant un premier groupe de 6 pays ayant pour mortalité (18,19,20,20,21,22) et un second groupe de 10 pays ayant pour mortalité (13,14,15,15,16,17). Les moyennes des deux groupes sont toujours respectivement de 20 et de 15, mais on “sent” que les écarts sont plus pertinents car ils portent sur deux groupes plus larges et surtout plus homogènes. Dit autrement, il est peu probable que les différences observées dans la seconde expérience soient le fruit du hasard alors que, dans le premier cas, on pouvait imaginer que les pays des deux groupes n’étaient pas fondamentalement différents.

# Comparaison de la distribution de deux échantillons
  VARQUANT<-c(18,19,20,20,21,22,13,14,15,15,16,17)
  VARQUALI<-c("A","A","A","A","A","A","B","B","B","B","B","B")
  par(mfrow=c(1,1))
  boxplot(VARQUANT~ VARQUALI)

Le test T de student dont vous trouverez la formule dans tous les bons manuels va permettre de transformer l’intuition précédente en mesure objective et en donner une mesure fondamentale appelée p-value ou seuil de rejet de l’hypothèse d’indépendance entre deux caractères. Là encore, on peut donner une idée intuitive de ce qu’est la p-value en la présentant comme une mesure de la significativité de la relation entre deux caractères. De façon plus rigoureuse, la p-value mesure la probabilité que la relation observée entre deux variables (ici, la différence de moyenne) soit l’effet du hasard. La p-value varie entre 0 et 1, 0 signifiant que la significativité de la relation est certaine (la différence de moyenne ne peut pas être l’effet du hasard) et 1 signifiant qu’il n’existe aucune relation entre les deux variobles (les deux moyennes sont égales).

Examinons brièvement les valeurs de p-value sur nos deux exemples :

# Test paramétrique d'égalité des moyennes de Student
  VARQUANT<-c(10,20,30,5,15,25)
  VARQUALI<-c("A","A","A","B","B","B")
   t.test(VARQUANT ~ VARQUALI) 

      Welch Two Sample t-test

  data:  VARQUANT by VARQUALI
  t = 0.61237, df = 4, p-value = 0.5734
  alternative hypothesis: true difference in means is not equal to 0
  95 percent confidence interval:
   -17.66958  27.66958
  sample estimates:
  mean in group A mean in group B 
               20              15 

Vous pouvez lire dans le résultat que la p-value est de 0.5734 ce qui signifie que l’on aurait 57% de chances de se tromper si l’on affirmait que la différence de moyenne entre les deux échantillons A et B est l’effet du hasard. Voyons ce qu’il en est pour notre second exemple

# Test paramétrique d'égalité des moyennes de Student
  VARQUANT<-c(18,19,20,20,21,22,13,14,15,15,16,17)
  VARQUALI<-c("A","A","A","A","A","A","B","B","B","B","B","B")
  t.test(VARQUANT ~ VARQUALI) 

      Welch Two Sample t-test

  data:  VARQUANT by VARQUALI
  t = 6.1237, df = 10, p-value = 0.0001121
  alternative hypothesis: true difference in means is not equal to 0
  95 percent confidence interval:
   3.180732 6.819268
  sample estimates:
  mean in group A mean in group B 
               20              15 

Cette fois-ci, vous pouvez voir que la valuer du t de student est beaucoup plus forte (6.13) ce qui se traduit par une p-value très proche de zéro (0.00011). La probabilité que les différences entre les deux échantillons soit l’effet du hasard est donc minime et on peut conclure avec une quasi-certiude que A et B correpondent à des groupes de pays ayant des mortalité significativement différentes.

Essayons maintenant d’appliquer ce même test de Student à notre exemple empirique de la mortalité infantile des pays socialistes et capitalistes.

# Test paramétrique d'égalité des moyennes de Student (inadapté)
   t.test(euro$TMI ~ euro$BLOC) 

      Welch Two Sample t-test

  data:  euro$TMI by euro$BLOC
  t = -3.1946, df = 7.3701, p-value = 0.01417
  alternative hypothesis: true difference in means is not equal to 0
  95 percent confidence interval:
   -21.06338  -3.24912
  sample estimates:
  mean in group Cap mean in group Soc 
            9.06875          21.22500 

Cette fois-ci le t de Student est égal à -3.19 et la p-value est égale à 0.0147, ce qui signifie que nous avons environ 1.5% de chances que les différences observées entre pays socialistes et capitalistes soient le simple effet du hasard et ne reflètent pas une véritable opposition. Du coup, que faut-il conclure ?

En sciences sociales, il est de tradition de considérer comme significatives les relations comportant un risque d’erreur inférieur à 5% soit 0.05. Puisque notre p-value est inférieure à 0.05, nous pouvons donc nous appuyer sur cette tradition et écrire que “l’on peut affirmer avec un risque d’erreur inférieur à 5% que les pays socialistes et capitalistes présentaient des différences significatives de mortalité infantile en 1988”.

Ou bien écrire plus simplement “En 1988 en Europe, la mortalité infantile moyenne des pays socialistes (21.2 P.1000) était significativement plus forte que celles des pays capitalistes (9.1 p.1000)” et ajouter dans une note en bas de page la remarque qui fait chic p-value=0.0147.

Ne concluez cependant pas trop vite que vous maitrisez désormais toutes les subtilités du test statistique d’égalité des moyennes. Si notre raisonnement précédent est approximativement juste, il serait sans doute critiqué par un vrai statisticien qui noterait immédiatement que le test t de student n’est pas le plus adapté à notre problème car (1) les distributions de mortalité ne sont pas gaussiennes et comportent des valeurs exceptionnelles comme l’Albanie pour les pays socialites ou le Portugal pour les pays capitalistes et (2) la variance des deux échantillons (leur hétérogénéité interne) n’est pas égale ce qui peut fausser l’usage du test de Student. Dans le cas présent, un bon statisticien préférerait sans doute utiliser un test non paramétrique, fondé sur la comparaison de l’ordre des observations et non pas leur valeur.

Soit par exemple le test de Wilcoxon, qu’il est très facile de mettre en oeuvre sour R en changeant à peine notre programme.

# Test non paramétrique d'égalité des médianes de Wilcoxon (correct)
   wilcox.test(euro$TMI~euro$BLOC) 

      Wilcoxon rank sum exact test

  data:  euro$TMI by euro$BLOC
  W = 8, p-value = 0.0001822
  alternative hypothesis: true location shift is not equal to 0

LA p-value est désormais égale à 0.0001 et on peut conclure avec certitude à l’existence de différence entre nos deux échantillons, ce qui n’était pas le cas avec le test de student, moins adapté au cas étudié.

2.4.1 Calcul des coefficients de corrélations

La façon la plus simple de calculer un coefficient de corrélation est d’appliquerla fonction cor à un tableau ne contenant que des variables quantitatives

# Coefficients de corrélation de pearson
  tab<-euro[,c("TMI","PNB","URB","FEC","POP")]
  cor(tab)

           TMI         PNB        URB        FEC         POP
  TMI  1.0000000 -0.67645013 -0.6523693  0.8141417 -0.13586790
  PNB -0.6764501  1.00000000  0.4929530 -0.6045034  0.01784716
  URB -0.6523693  0.49295300  1.0000000 -0.5436824  0.39790430
  FEC  0.8141417 -0.60450337 -0.5436824  1.0000000 -0.20123367
  POP -0.1358679  0.01784716  0.3979043 -0.2012337  1.00000000

R calcule par débaut le coefficient de corrélation linéaire de Bravais-Pearson qui varie entre -1 (corrélation négative parfaite) et +1 (corrélation positive parfaite). On peut ainsi voir que le taux de mortalité infantile (TMI) a une très forte corrélation positive avec l’indice conjonctuel de fécondité (+0.81) et une forte corrélation négative avec le PNB/hab (-0.68) ou le taux d’urbanisation (-0.65). La corrélationavec la population du pays est négative mais faible (-0.14).

2.4.2 Test de significativité d’un coefficient de corrélation

Cette procédure n’est cependant pas totalement convaincante d’un point de vue statistique car elle n’affiche pas le degré de significativité de la liaison statistique, la fameuse p-value dont nous avons discuté dans l’exemple précédent. Cette significativité dépend en effet à la fois (1) de la valeur absolue du coefficient de corrélation (plus on se rapproche de +1 ou -1,plus la relation est significative) mais aussi (2) de la taille de l’échantillon utilisé pour calculer cette corrélation, appelée nombre de degrés de liberté. Lorsqu’on calcule un coefficient de corrélation, le nombre de dégré de liberté est égal à (N-2) c’est-àdire au nombre d’observations (N) moins le nombre de paramètres utilisés pour faire le calcul (soit 2, puisqu’on a besoin de la moyenne de X et de celle de Y).

Pour savoir si une relation est significative, il faut procéder à un test statistique à l’aide de la fonction cor.test. On va par exemple tester la significativité des relations entre TMI et les quatre autres variables.

# Test du coefficient de corrélation de Peason
  
  cor.test(tab$TMI,tab$PNB, type = "pearson")

      Pearson's product-moment correlation

  data:  tab$TMI and tab$PNB
  t = -4.3081, df = 22, p-value = 0.0002843
  alternative hypothesis: true correlation is not equal to 0
  95 percent confidence interval:
   -0.8483508 -0.3755262
  sample estimates:
         cor 
  -0.6764501 

cor.test(tab$TMI,tab$URB, type = "pearson")

      Pearson's product-moment correlation

  data:  tab$TMI and tab$URB
  t = -4.0373, df = 22, p-value = 0.0005507
  alternative hypothesis: true correlation is not equal to 0
  95 percent confidence interval:
   -0.835811 -0.337894
  sample estimates:
         cor 
  -0.6523693 

cor.test(tab$TMI,tab$FEC, type = "pearson")

      Pearson's product-moment correlation

  data:  tab$TMI and tab$FEC
  t = 6.5763, df = 22, p-value = 1.296e-06
  alternative hypothesis: true correlation is not equal to 0
  95 percent confidence interval:
   0.6116118 0.9165298
  sample estimates:
        cor 
  0.8141417 

cor.test(tab$TMI,tab$POP, type = "pearson")

      Pearson's product-moment correlation

  data:  tab$TMI and tab$POP
  t = -0.64324, df = 22, p-value = 0.5267
  alternative hypothesis: true correlation is not equal to 0
  95 percent confidence interval:
   -0.511244  0.283042
  sample estimates:
         cor 
  -0.1358679 

On voit ainsi apparaître deux tableaux, l’un avec les coefficients de corrélation des variables et l’autre avec les p-value associées à chacun de ces coefficients. Il apparaît alors que les corrélations du taux de mortalité infantile sont très significatives (<0.001) avec le PNB/hab, le taux d’urbanisation et l’indice conjoncturel de fécondité. Mais que la corrélation entre mortalité infantile et population est non-significative.

Ajoutons maintenant une autre précaution statistique consistant à calculer un autre coefficient de corrélation, le coefficient de corrélation de rang aussi appelé coefficient de Spearman. L’intérêt de ce coefficient est tout d’abord de limiter le jeu des valeurs exceptionelles (outlier) qui sont liées à la position excentré d’une seule observation. Mais c’est aussi plus généralement de mettre en évidence des relations croissantes ou décroissantes de forme non-linéaire entre deux variables. En présence de ces deux types de singularité, le coefficient de Spearman affiche souvent des valeurs très divergentes de celui de Pearson. Et la découverte d’une discordance entre les deux coefficients est un indice à ne pas négliger.

# Coefficient de corrélation de Spearman
  cor(tab, method = "spearman")

           TMI        PNB        URB         FEC         POP
  TMI  1.0000000 -0.9018927 -0.5690669  0.50791428  0.16782609
  PNB -0.9018927  1.0000000  0.5019591 -0.64715761 -0.18533830
  URB -0.5690669  0.5019591  1.0000000 -0.43758300  0.35196875
  FEC  0.5079143 -0.6471576 -0.4375830  1.00000000 -0.07105537
  POP  0.1678261 -0.1853383  0.3519688 -0.07105537  1.00000000

# Test du cofficient de Spearman
  cor.test(tab$TMI,tab$PNB, method = "spearman")

      Spearman's rank correlation rho

  data:  tab$TMI and tab$PNB
  S = 4374.4, p-value = 1.763e-09
  alternative hypothesis: true rho is not equal to 0
  sample estimates:
         rho 
  -0.9018927 

cor.test(tab$TMI,tab$URB, method = "spearman")

      Spearman's rank correlation rho

  data:  tab$TMI and tab$URB
  S = 3608.9, p-value = 0.003707
  alternative hypothesis: true rho is not equal to 0
  sample estimates:
         rho 
  -0.5690669 

cor.test(tab$TMI,tab$FEC, method = "spearman")

      Spearman's rank correlation rho

  data:  tab$TMI and tab$FEC
  S = 1131.8, p-value = 0.01128
  alternative hypothesis: true rho is not equal to 0
  sample estimates:
        rho 
  0.5079143 

cor.test(tab$TMI,tab$POP, method = "spearman")

      Spearman's rank correlation rho

  data:  tab$TMI and tab$POP
  S = 1914, p-value = 0.4314
  alternative hypothesis: true rho is not equal to 0
  sample estimates:
        rho 
  0.1678261 

Dans l’exemple (ancien mais pédagogique) qui a été choisi, il y a clairement une singularité dans la relation entre TMI et PNB puisque la corrélation qui était de (+0.68) avec le coefficient de Pearson grimpe à (+0.90) avec celui de Spearman. Dans le même temps, la relation entre TMI et FEC chute de (+0.81) à (+0.51) ce qui indique une autre singularité. Il faudra donc examiner soigneusement les nuages de points pour comprendre d’où viennent ces divergences, ce que l’on va faire au cours de l’étape suivante.

2.4.3 Visualisation de la forme de la relation entre deux variables quantitatives

On peut visualiser une relations entre deux variables X et Y à l’aide de l’instruction générique plot(). Soit par exemple la relation entre urbanisation et mortalité infantile :

plot(euro$URB,euro$TMI)

On peut également ranger plusieurs graphiques en une seule figure en utilisant l’instruction générale par qui définit le format des graphiques, associée au paramètre mfrow qui indique combien de lignes et de colonnes on souhaite générer. Chacune des case de ce tableau servira successivement à afficher les graphiques demandés par l’instruction plot. Ici, nous déclarons que nous voulons ordonner les graphiques sur 2 lignes et 2 colonnes puis nous executons 4 fois de suite l’instruction plot.

# Un seul graphique de synthèse
  par(mfrow= c(2,2)) 
  plot(euro$PNB,euro$TMI)
  plot(euro$URB,euro$TMI)
  plot(euro$FEC,euro$TMI)
  plot(euro$POP,euro$TMI)

Ces figures permettent (avec un peu d’habitude …) d’expliquer les divergences observées entre les coefficients de corrélation de Pearson et de Spearman.

par(mfrow= c(1,2)) 
  
  X<-euro$PNB
  Y<-euro$TMI
  
  corlin<-round(cor(X,Y),3)
  reglin<-lm(Y~X)
  plot(X,Y, main=paste("Coeff Pearson = ",corlin))
  abline(reglin,col="red")
  
  rankX<-rank(X)
  rankY<-rank(Y)
  corrank<-round(cor(rankX,rankY),3)
  regrank<-lm(rankY~rankX)
  plot(rankX,rankY, main=paste("Coeff Spearman = ",corrank))
  abline(regrank,col="red")

• la relation entre TMI et PNB est clairement non linéaire et forme un arc de parabole plutôt qu’une droite. C’est la raison pour laquelle le coefficient de Spearman était beaucoup plus élevé que celui de Pearson. En effet, lorsque l’on transforme les variables en rang, la courbure disparaît.

par(mfrow= c(1,2)) 
  
  X<-euro$FEC
  Y<-euro$TMI
  
  corlin<-round(cor(X,Y),3)
  reglin<-lm(Y~X)
  plot(X,Y, main=paste("Coeff Pearson = ",corlin))
  abline(reglin,col="red")
  
  rankX<-rank(X)
  rankY<-rank(Y)
  corrank<-round(cor(rankX,rankY),3)
  regrank<-lm(rankY~rankX)
  plot(rankX,rankY, main=paste("Coeff Spearman = ",corrank))
  abline(regrank,col="red")

plot(euro$URB,euro$TMI)

On peut également ranger plusieurs graphiques en une seule figure en utilisant l’instruction générale par qui définit le format des graphiques, associée au paramètre mfrow qui indique combien de lignes et de colonnes on souhaite générer. Chacune des case de ce tableau servira successivement à afficher les graphiques demandés par l’instruction plot. Ici, nous déclarons que nous voulons ordonner les graphiques sur 2 lignes et 2 colonnes puis nous executons 4 fois de suite l’instruction plot.

# Un seul graphique de synthèse
  par(mfrow= c(2,2)) 
  plot(euro$PNB,euro$TMI)
  plot(euro$URB,euro$TMI)
  plot(euro$FEC,euro$TMI)
  plot(euro$POP,euro$TMI)

Ces figures permettent (avec un peu d’habitude …) d’expliquer les divergences observées entre les coefficients de corrélation de Pearson et de Spearman.

par(mfrow= c(1,2)) 
  
  X<-euro$PNB
  Y<-euro$TMI
  
  corlin<-round(cor(X,Y),3)
  reglin<-lm(Y~X)
  plot(X,Y, main=paste("Coeff Pearson = ",corlin))
  abline(reglin,col="red")
  
  rankX<-rank(X)
  rankY<-rank(Y)
  corrank<-round(cor(rankX,rankY),3)
  regrank<-lm(rankY~rankX)
  plot(rankX,rankY, main=paste("Coeff Spearman = ",corrank))
  abline(regrank,col="red")

• la relation entre TMI et PNB est clairement non linéaire et forme un arc de parabole plutôt qu’une droite. C’est la raison pour laquelle le coefficient de Spearman était beaucoup plus élevé que celui de Pearson. En effet, lorsque l’on transforme les variables en rang, la courbure disparaît.

par(mfrow= c(1,2)) 
  
  X<-euro$FEC
  Y<-euro$TMI
  
  corlin<-round(cor(X,Y),3)
  reglin<-lm(Y~X)
  plot(X,Y, main=paste("Coeff Pearson = ",corlin))
  abline(reglin,col="red")
  
  rankX<-rank(X)
  rankY<-rank(Y)
  corrank<-round(cor(rankX,rankY),3)
  regrank<-lm(rankY~rankX)
  plot(rankX,rankY, main=paste("Coeff Spearman = ",corrank))
  abline(regrank,col="red")

• la relation entre TMI et FEC est quant à elle influencée par une valeur exceptionnelle, l’Albanie, qui affiche une fécondité et une mortalité infantile excpetionnnelles. Le point Albanie est situé très en dehors du nuage des autres points et il fait fortement augmenter le coefficient de Pearson. Mais si l’on passe en rang, son effet disparaît et le nuage de point apparaît beaucoup moins linéaire, ce qui explique la valeur plus faible du coefficient de Spearman.

Le point important à retenir (en l’absence de la lecture d’un manuel de statistique plus précis) est qu’on ne doit pas tenter de faire passer une droite dans un nuage de points si celui-ci ne forme pas un ensemble de points régulièrement alignés. Il faut éviter d’utiliser la régression linéaire sur un nuage de points ayant une courbure ou sur un nuage de points comportant des valeurs exceptionnelles. On va voir dans la section suivante quelques “trucs” permettant de corriger la forme d’un nuage de point inadapté à une analyse de régression linéaire (ce qui ne dispense pas de lire un manuel de statistique expliquant ces trucs à l’aide des concepts plus précis d’autocorrélation des résidus, d’hétéroscédasticité, etc…)

2.4.4 Un premier modèle de régression linéaire (quick & dirty …)

Le calcul d’une régression linéaire se fait très facilement avec l’instruction lm qui est l’acronyme de linear model. La particularité de R par rapport à d’autres logiciels est de stocker les résultats d’une régression (et plus généralement d’un modèle statistique quelconque) dans un objet dont on choisit le nom (ici : MonModele). C’est objet est uen sorte de grand sac sans fonds contenant des dizaines de résultats que l’on peut sortir un par un en fonction de ses besoins et surtout de ses connaissances en statistiques… Les débutants se limiteront à demander de sortir du sac à malice un résumé (summary) de la régression qu’ils viennent d’opérer :

MonModele <- lm(euro$TMI~euro$PNB)
  summary(MonModele)

  Call:
  lm(formula = euro$TMI ~ euro$PNB)

  Residuals:
     Min     1Q Median     3Q    Max 
  -7.921 -3.222 -1.440  1.064 22.344 

  Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
  (Intercept) 21.3403638  2.3136443   9.224 5.14e-09 ***
  euro$PNB    -0.0011403  0.0002647  -4.308 0.000284 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

  Residual standard error: 6.411 on 22 degrees of freedom
  Multiple R-squared:  0.4576,    Adjusted R-squared:  0.4329 
  F-statistic: 18.56 on 1 and 22 DF,  p-value: 0.0002843

On trouve dans ce résumé l’essentiel à savoir l’équation de la droite de régression linéaire dans la partie coefficient (TMI = -0.00114 * PNB + 21.34), le pouvoir explicatif du modèle qui est le carré du coefficient de corrélation appelé en anglais R-Squared (0.457 soit 46%) et la significativité générale du modèle qui est fournie par la p-value (0.00028). Cette dernière est en l’occurence conforme au coefficient de corrélation entre TMI et PNB car nous faisons une régression simple. On verra ultérieurement le cas des régressions multiples.

La visualisation de cette droite de régression peut se faire très simplement en traçant le nuage de point avec plot et en ajoutant la droite en utilisant l’instruction abline appliquée au sac à malice où se trouvent les résultats de notre modèle. Même si nous ignorons le contenu détaillé de ce sac, nous savons que R va y trouver les éléments utiles à savoir les coefficients de la droite de régression que nous avons déjà vu dans le résumé. On va faire un petit effort d’habillage en ajoutant à la figure un titre principal (main) et des titres d’axes (xlab,ylab).

par(mfrow=c(1,1))
  plot(euro$PNB,euro$TMI,
       main=" Richesse et mortalité infantile en Europe vers 1988",
       ylab="PNB en $ par habitant (1988)",
       xlab="Taux de mortalité infantile en p. 1000 (1988)") 
  abline(MonModele, col="red")

par(mfrow=c(1,1))
  plot(euro$PNB,euro$TMI,
       main=" Richesse et mortalité infantile en Europe vers 1988",
       ylab="PNB en $ par habitant (1988)",
       xlab="Taux de mortalité infantile en p. 1000 (1988)") 
  abline(MonModele, col="red")

Nous avons réussi à faire notre premier modèle de régression sous R ! Mais la vérité oblige à dire qu’il n’est pas satisfaisant puisque nous avons ajusté une droite dans un nuage de point qui est manifestement de forme non linéaire. Essayons de faire mieux, en montrant que la question statistique peut recouper des question plus théoriques de géographie économique et de géographie politique.

2.4.5 Un deuxième modèle plus satisfaisant sur le plan statistique et théorique

La courbure de la relation entre PNB/habitant et taux de mortalité infantile est certes un problème statistique. Mais elle soulève aussi une question théorique plus intéressante : Est-il raisonnable et logique de penser que plus le PNB/hab. augmente, plus la mortalité infantiel décroît de façon linéaire ?

Notre modèle précédent (Y=aX+b) est en partie absurde puisqu’il nous dit que la mortalité infantile minimale est de 21.34°/°° dans un pays ayant un PNB nul (on sait qu’il existe des mortalité infantile plus forte) et surtout il suppose que chaque fois que le PNB augmente de 1000 $, la mortalité infantile baisse de 1.1 °/°°. Du coup, ce modèle prévoit une absence totale de mortalité infantile pour les pays ayant un revenu supérieur à 15000 $/hab. Et il annonce froidement des mortalités infantile négatives au delà de ce niveau de richesse ! Bref, il y a une inconsistance logique dans un tel modèle, indépendamment du problème statistique de non linéarité.

Or, ce que nous disent les théoriciens du développement comme A. Sen est l’existence de rendements décroissant dans les progrès au fur et à mesure que la richesse augmente. En d’autres termes, la mortalité infantile diminue rapidement quant on passe de l’extrême apuvreté à la pauvreté. Mais ensuite elle baisse de moins en moins vite et son effet devient négligeable au delà d’un seuil de richesse. Le modèle correspondant à ces hypothèses théoriques est donc un modèle de décroissance non linéaire.

Il peut s’agir d’une loi exponentielle négative qui peut s’écrire Y = b.exp(-aX), ce que l’on peut transformer en équation linéaire en le réécrivant sous la forme log(Y) = log(b)+aX. Ou bien d’une loi puissance négative qui peut s’écrire Y = b.X^a et se transforme en équation linéaire log(y) = a.log(X)+log(b). Si l’on n’a pas de préférence théorique pour l’une ou l’autre de ces lois, on peut simplement examiner celle qui donne le meilleur ajustement linéaire après transformation logarithmique :

par(mfrow= c(1,2)) 
  
  X<-euro$PNB
  Y<-euro$TMI
  
  corexp<-round(cor(X,log(Y)),3)
  regexp<-lm(log(Y)~X)
  plot(X,log(Y), main=paste("Modèle exponentiel : r = ",corexp))
  abline(regexp,col="red")
  
  corpui<-round(cor(log(X),log(Y)),3)
  regpui<-lm(log(Y)~log(X))
  plot(log(X),log(Y), main=paste("Modèle puissance : r = ",corpui))
  abline(regpui,col="red")

par(mfrow= c(1,2)) 
  
  X<-euro$PNB
  Y<-euro$TMI
  
  corexp<-round(cor(X,log(Y)),3)
  regexp<-lm(log(Y)~X)
  plot(X,log(Y), main=paste("Modèle exponentiel : r = ",corexp))
  abline(regexp,col="red")
  
  corpui<-round(cor(log(X),log(Y)),3)
  regpui<-lm(log(Y)~log(X))
  plot(log(X),log(Y), main=paste("Modèle puissance : r = ",corpui))
  abline(regpui,col="red")

La comparaison des coefficients de corrélation mais aussi l’examen de la forme des nuages de points plaide clairement en faveur du choix du modèle puissance. Il demeure en effet une nette courbure du nuage de points dans le graphique du modèle exponentiel et l’ajustement obtenu est plus faible. Le modèle puissance affiche un alignement presque parfait des points avec une corrélation très élevée. Donc, si l’on ne dispose pas d’arguments théoriques ou empiriques particuliers en faveur du modèle exponentiel, on retiendra le modèle puissance pour modéliser la relation entre mortalité infantile et PNB/habitant des pays européens en 1988.

2.4.6 Un troisième modèle : la revanche des pays socialistes …

Comme nous savons désormais (1) que la mortalité infantile varie fortement entre les pays socialistes et capitalistes et (2) que la mortalité infantile dépend de la richesse par habitant, on peut proposer un dernier modèle qui va faire intervenir simultanément les variables. Il existe de nombreuses variantes pour construire un tel modèle que vous apprendrez dans les cours sur la régression multiple ou l’analyse multi-niveaux. On se limitera ici à une approche plus intuitive qui ne présuppose pas encore de connaître des techniques statistiques poussées.

Examinons tout d’abord la relation TMI*PNB dans les pays socialistes :

par(mfrow=c(1,1))
  summary(euro$PNB)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      600    2275    6850    7208   10575   17800 

eurosoc<-euro[euro$BLOC=="Soc",]
  MonModele2soc <- lm(eurosoc$TMI~eurosoc$PNB)
  summary(MonModele2soc)

  Call:
  lm(formula = eurosoc$TMI ~ eurosoc$PNB)

  Residuals:
      Min      1Q  Median      3Q     Max 
  -7.9609 -4.0541 -0.7209  3.5026  7.9557 

  Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
  (Intercept) 40.437252   5.552599   7.283 0.000341 ***
  eurosoc$PNB -0.008988   0.002383  -3.772 0.009271 ** 
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

  Residual standard error: 6.25 on 6 degrees of freedom
  Multiple R-squared:  0.7033,    Adjusted R-squared:  0.6539 
  F-statistic: 14.22 on 1 and 6 DF,  p-value: 0.009271

plot(eurosoc$PNB,eurosoc$TMI,
                             xlim=c(0,6000),
                              ylim=c(0,50),
                             main=" Pays socialistes",
                             ylab="PNB en $ par habitant (1988)",
                             xlab="Taux de mortalité infantile en p. 1000 (1988)") 
  abline(MonModele2soc, col="red")

plot(eurosoc$PNB,eurosoc$TMI,
                             xlim=c(0,6000),
                              ylim=c(0,50),
                             main=" Pays socialistes",
                             ylab="PNB en $ par habitant (1988)",
                             xlab="Taux de mortalité infantile en p. 1000 (1988)") 
  abline(MonModele2soc, col="red")

La relation semble bien linéaire avec une valeur initiale forte de mortalité infantile (coefficient b=40.4) et une décroissance rapide de la mortalité infantile lorsque le PNB augmente (coefficient a = -0.009). La mortalité infantile baisse de 9 points chaque fois que le PNB augmente de 1000 $. On peut estimer qu’au delà d’un PNB de 5000 $ environ, la mortalité infantile n’existera plus dans les pays socialistes.

Examinons maintenant séparément la situation dans les pays capitalistes

par(mfrow=c(1,1))
  eurocap<-euro [euro$BLOC=="Cap",]
  MonModele2cap <- lm(eurocap$TMI~eurocap$PNB)
  summary(MonModele2cap)

  Call:
  lm(formula = eurocap$TMI ~ eurocap$PNB)

  Residuals:
      Min      1Q  Median      3Q     Max 
  -2.3212 -1.3747  0.4413  0.7973  3.5595 

  Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
  (Intercept) 13.1654634  1.0987927  11.982 9.54e-09 ***
  eurocap$PNB -0.0004204  0.0001039  -4.048   0.0012 ** 
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

  Residual standard error: 1.711 on 14 degrees of freedom
  Multiple R-squared:  0.5392,    Adjusted R-squared:  0.5063 
  F-statistic: 16.38 on 1 and 14 DF,  p-value: 0.001199

plot(eurocap$PNB,eurocap$TMI,
      xlim=c(0,20000),
      ylim=c(0,20),
       main=" Pays capitalistes",
       ylab="PNB en $ par habitant (1988)",
       xlab="Taux de mortalité infantile en p. 1000 (1988)") 
  abline(MonModele2cap, col="blue")

La relation semble à nouveauà peu près linéaire. La valeur initiale de mortalité infantile est plus faible que dans les pays socialistes (coefficient a = 13.1) mais ensuite la décroissance est lente lorsque le PNB augmente (coefficient a = -0.0004). La mortalité infantile ne baisse que de 0.4 points lorsque le PNB augmente de 1000 $. Il faudrait donc un PNB d’environ 30000 $ par habitant pour que la mortalité infantile disparaissent dans les pays capitalistes, alors qu’il suffisant de 5000 $ par habitant dans les pays socialistes !

On va résumer ce dernier modèle par une seule figure sur laquelle on affichera simultanément les deux droites de régression.

par(mfrow=c(1,1))
  
  plot(euro$PNB,euro$TMI,
      xlim=c(0,20000),
      ylim=c(0,50),
       main="Influence du niveau de vie et du syème politique sur la mortalité infantile",
       xlab="PNB en $ par habitant (1988)",
       ylab="Taux de mortalité infantile en p. 1000 (1988)") 
  
  points(eurosoc$PNB,eurosoc$TMI,col="red")
  text(eurosoc$PNB,eurosoc$TMI,eurosoc$PAYS,col="red",cex=0.5,pos=4)
  abline(MonModele2soc,col="red")
  
  points(eurocap$PNB,eurocap$TMI,col="blue")
  text(eurocap$PNB,eurocap$TMI,eurocap$PAYS,col="blue",cex=0.5,pos=4)
  abline(MonModele2cap,col="blue")

plot(eurocap$PNB,eurocap$TMI,
      xlim=c(0,20000),
      ylim=c(0,20),
       main=" Pays capitalistes",
       ylab="PNB en $ par habitant (1988)",
       xlab="Taux de mortalité infantile en p. 1000 (1988)") 
  abline(MonModele2cap, col="blue")

La relation semble à nouveauà peu près linéaire. La valeur initiale de mortalité infantile est plus faible que dans les pays socialistes (coefficient a = 13.1) mais ensuite la décroissance est lente lorsque le PNB augmente (coefficient a = -0.0004). La mortalité infantile ne baisse que de 0.4 points lorsque le PNB augmente de 1000 $. Il faudrait donc un PNB d’environ 30000 $ par habitant pour que la mortalité infantile disparaissent dans les pays capitalistes, alors qu’il suffisant de 5000 $ par habitant dans les pays socialistes !

On va résumer ce dernier modèle par une seule figure sur laquelle on affichera simultanément les deux droites de régression.

par(mfrow=c(1,1))
  
  plot(euro$PNB,euro$TMI,
      xlim=c(0,20000),
      ylim=c(0,50),
       main="Influence du niveau de vie et du syème politique sur la mortalité infantile",
       xlab="PNB en $ par habitant (1988)",
       ylab="Taux de mortalité infantile en p. 1000 (1988)") 
  
  points(eurosoc$PNB,eurosoc$TMI,col="red")
  text(eurosoc$PNB,eurosoc$TMI,eurosoc$PAYS,col="red",cex=0.5,pos=4)
  abline(MonModele2soc,col="red")
  
  points(eurocap$PNB,eurocap$TMI,col="blue")
  text(eurocap$PNB,eurocap$TMI,eurocap$PAYS,col="blue",cex=0.5,pos=4)
  abline(MonModele2cap,col="blue")

Selon ce modèle, on pouvait penser en 1988 qu’un pays capitaliste comme le Portugal devrait très fortementaccroître son PNB par habitant pour voir baisser sa mortalité infantile. Inversement, la Pologne qui était un pays socialiste aurait du voir rapidement sa situation s’améliorer, pour peu qu’elle accroisse légèrement son PNB/habitant. Mais on ne saura jamais si le modèle aurait été vérifié puisque le bloc socialiste s’est dissous en 1989 …

3.1.1 Définir un espace de travail

On vérifie l’emplacement où l’on va trouver les données

# (1.2) Choix du rÈpertoire contenant les donnÈes
  monchemin<-"data/yemessoa"
  list.files(monchemin)

[1] "yemessoa2005.txt" "yemessoa2005.xls"

Nous allons maintenant essayer de charger le fichier de données

3.2.1 Importation d’un tableau

# --------------------------------------------------------
  # (2) IMPORTATION ET MISE EN FORME D'UN TABLEAU DE DONNEES
  # --------------------------------------------------------
  # (2.1) Importation d'un fichier .txt
  yem<-read.csv("data/yemessoa/yemessoa2005.txt", 
                header=TRUE, 
                sep="\t",
                dec=",",
                stringsAsFactors = T)

On affiche les premières lignes (head) ou les dernières lignes (tail) du tableau pour voir si le chargement s’est bien effectué :

head(yem)

  CODE CON2 CON3 DIST               QUAR          VIL      CLAN  SEXE      SELF
  1  237  Oui    2    3       Elig Messobo  Yemessoa II   B\xe9ti Homme (1) Petit
  2   98  Non    0   20 Nkol Ngu\xe9gu\xe9   Yemessoa I B\xe9loua Homme (1) Petit
  3  157  Non    0   30       Biling Bitom   Yemessoa I B\xe9loua Homme (1) Petit
  4  325  Oui    1   10       Elig Messobo  Yemessoa II   B\xe9ti Homme (2) Moyen
  5  299  Oui    1   10       Nkol Ambenbe Yemessoa III   B\xe9ti Homme (3) Grand
     BUD TABLE SALON AGE             ETU5                       ETU2
  1 3950   Non   Non  19   ETU3 : CEP-CAP ETU2 : Sup\xe9rieur au CEP
  2 3700   Oui   Non  20   ETU3 : CEP-CAP ETU2 : Sup\xe9rieur au CEP
  3 1250   Non   Non  21 ETU5 : BEPC et + ETU2 : Sup\xe9rieur au CEP
  4 2200   Oui   Non  21   ETU3 : CEP-CAP ETU2 : Sup\xe9rieur au CEP
  5 4450   Oui   Non  21   ETU3 : CEP-CAP ETU2 : Sup\xe9rieur au CEP
   [ reached 'max' / getOption("max.print") -- omitted 1 rows ]

tail(yem)

    CODE CON2 CON3 DIST          QUAR          VIL    CLAN  SEXE      SELF  BUD
  294  128  Oui    2    5  Biling Bitom   Yemessoa I B\xe9ti Homme (3) Grand 3200
  295   71  Non    0   35  Nkol Ambenbe Yemessoa III B\xe9ti Homme (3) Grand 2850
  296  176  Oui    2   15  Biling Bitom   Yemessoa I B\xe9ti Femme (3) Grand 2675
  297  178  Oui    2   10  Biling Bitom   Yemessoa I B\xe9ti Homme (3) Grand 3750
  298  243  Oui    2   25 Mebang Mengoe   Yemessoa I B\xe9ti Homme (3) Grand 8900
      TABLE SALON AGE                      ETU5                      ETU2
  294   Oui   Non  81               ETU1: Aucun ETU1: inf\xe9rieur au CEP
  295   Oui   Non  83 ETU2 :  Primaire sans CEP ETU1: inf\xe9rieur au CEP
  296   Oui   Oui  87               ETU1: Aucun ETU1: inf\xe9rieur au CEP
  297   Oui   Oui  88 ETU2 :  Primaire sans CEP ETU1: inf\xe9rieur au CEP
  298   Oui   Oui  90 ETU2 :  Primaire sans CEP ETU1: inf\xe9rieur au CEP
=======
  
head(yem)
  
  CODE CON2 CON3 DIST         QUAR          VIL   CLAN  SEXE      SELF  BUD
  1  237  Oui    2    3 Elig Messobo  Yemessoa II   Béti Homme (1) Petit 3950
  2   98  Non    0   20 Nkol Nguégué   Yemessoa I Béloua Homme (1) Petit 3700
  3  157  Non    0   30 Biling Bitom   Yemessoa I Béloua Homme (1) Petit 1250
  4  325  Oui    1   10 Elig Messobo  Yemessoa II   Béti Homme (2) Moyen 2200
  5  299  Oui    1   10 Nkol Ambenbe Yemessoa III   Béti Homme (3) Grand 4450
    TABLE SALON AGE             ETU5                    ETU2
  1   Non   Non  19   ETU3 : CEP-CAP ETU2 : Supérieur au CEP
  2   Oui   Non  20   ETU3 : CEP-CAP ETU2 : Supérieur au CEP
  3   Non   Non  21 ETU5 : BEPC et + ETU2 : Supérieur au CEP
  4   Oui   Non  21   ETU3 : CEP-CAP ETU2 : Supérieur au CEP
  5   Oui   Non  21   ETU3 : CEP-CAP ETU2 : Supérieur au CEP
   [ reached 'max' / getOption("max.print") -- omitted 1 rows ]
  
tail(yem)
  
    CODE CON2 CON3 DIST          QUAR          VIL CLAN  SEXE      SELF  BUD
  294  128  Oui    2    5  Biling Bitom   Yemessoa I Béti Homme (3) Grand 3200
  295   71  Non    0   35  Nkol Ambenbe Yemessoa III Béti Homme (3) Grand 2850
  296  176  Oui    2   15  Biling Bitom   Yemessoa I Béti Femme (3) Grand 2675
  297  178  Oui    2   10  Biling Bitom   Yemessoa I Béti Homme (3) Grand 3750
  298  243  Oui    2   25 Mebang Mengoe   Yemessoa I Béti Homme (3) Grand 8900
      TABLE SALON AGE                      ETU5                   ETU2
  294   Oui   Non  81               ETU1: Aucun ETU1: inférieur au CEP
  295   Oui   Non  83 ETU2 :  Primaire sans CEP ETU1: inférieur au CEP
  296   Oui   Oui  87               ETU1: Aucun ETU1: inférieur au CEP
  297   Oui   Oui  88 ETU2 :  Primaire sans CEP ETU1: inférieur au CEP
  298   Oui   Oui  90 ETU2 :  Primaire sans CEP ETU1: inférieur au CEP
>>>>>>> d676e259c8be1e9f88c818dd0e1c0cc0a591abdd
   [ reached 'max' / getOption("max.print") -- omitted 1 rows ]